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Abstract 

We present results for the dependence of the residual mass of domain wall fermions (DWF) on 
the size of the fifth dimension and its relation to the density and localization properties of low-lying 
eigenvectors of the corresponding hermitian Wilson Dirac operator relevant to simulations of 2+1 
flavor domain wall QCD. Using the DBW2 and Iwasaki gauge actions, we generate ensembles of 
configurations with a 16^ x 32 space-time volume and an extent of 8 in the fifth dimension for the 
sea quarks. We demonstrate the existence of a regime where the degree of locality, the size of chiral 
symmetry breaking and the rate of topology change can be acceptable for inverse lattice spacings 
> 1.6 GeV. 

PACS numbers: 11.15.Ha, ll.30.Rd, 12.38.Aw, 12.38.-t 12.38.Gc 
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I. INTRODUCTION 



In this paper we present a study of three closely related quantities that are important 
for the successful application of the domain wall fermion (DWF) formulation to full QCD 
calculations with 2+1 flavors of light quarks. The first of these is the size of the explicit 
chiral symmetry breaking that occurs for domain wall fermions due to the finite extent of 
the lattice in the fifth dimension. This is usually characterized by the residual mass, mres, 
whose dependence on the gauge action, gauge coupling and the extent of the fifth dimension 
is discussed. 

The second topic covers the spectral properties of the four- dimensional hermitian Wilson 
Dirac operator evaluated at large, negative mass, — M5. All current formulations of lattice 
chiral symmetry make use of a negative-mass (hermitian) Wilson Dirac operator, Hw = 
75DVI/, and intellectually derive from the Kaplan approach [l], whether by taking an analytic 
limit as in the overlap operator P], or by directly simulating a finite fifth dimension 3] of 
size Lg. Auxiliary information on the nature of the spectrum of Hw is required to have 
confidence that a simulation at fixed lattice spacing lies in the universality class of QCD. 
This topic is closely connected to the properties of the Aoki phase of the related operator 
Dw [4]. If a chiral fermion calculation is attempted for a gauge action and a negative fermion 
mass that is too close to the Aoki phase, undesirable non-locality or enhanced explicit chiral 
symmetry breaking may result. Thus, we will examine both the density of the low-lying 
eigenmodes of Hw and their locality properties. 

Finally, we examine the ergodicity of the sampling of topological charge that results for 
these choices of gauge action and negative Wilson Dirac mass using the rational hybrid 



Monte Carlo method of Clark and Kennedy 



6|. Since a change in the topological charge 



of an evolving lattice configuration is expected to be accompanied by a zero mode of Hw, a 
choice of the gauge action which suppresses such zero modes will inhibit topology change, 
possibly rendering the Monte Carlo sampling non-ergodic in practice. 

Based on these studies, we conclude that lattice QCD calculations with 2+1 flavors of 
light quarks are possible with adequate locality, small explicit chiral symmetry breaking and 
adequate sampling of topological charge for inverse lattice spacings of 1.6 GeV or larger. 
In particular, the rapid decrease in the rate of change of topological charge for the DBW2 
action with decreasing lattice spacing favours the choice of Iwasaki action — a choice that has 
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been adopted for the large-scale ensemble generation of the RBC and UKQCD collaborations 
using the QCDOC computers [7|,l8|,l9| at the RIKEN-BNL Research Center and the University 
of Edinburgh. 

The structure of this paper is as follows. In Section [III we describe the Aoki phase of the 
Wilson Dirac operator in light of recent research and the connection to the spectrum and 
localization of eigenmodes of Hw- In Section UTI] we review the importance of the spectrum of 
Hw for lattice formulations that accurately realize chiral symmetry. The explicit breaking 
of chiral symmetry that results from finite Lg for domain wall fermions is discussed in 
Section [IV] and the expected dependence of the residual mass on Ls motivated. In Section IVl 
we describe the details of the gauge and DWF actions used in our simulations, while in 
Section |Vl] we present a microscopic study of the low-lying eigenvalues and eigenmodes of 
Hw on our ensembles. In Section IVIII we present a study of the dependence on Lg of the 
residual mass for valence quarks on each of the (fixed sea quark Lg) ensembles. The Lg 
dependence of these results are fitted to the form motivated in Section [IV] and conclusions 
about the low-lying spectrum of the transfer matrix in the fifth dimension are drawn and 
related to our results from Section IVII Finally, in Section IVIIIl we examine the dependence 
of the topological charge with evolving Monte Carlo time and demonstrate the correlation 
between reduced DWF chiral symmetry breaking and a suppression in the rate of change of 
topological charge. Appendix |X] summarizes the transfer matrix formalism of Furman and 
Shamir which is used extensively in the residual mass discussion of Section llVi 



II. AOKI PHASE DIAGRAM OF LATTICE QCD 

The properties of lattice chiral fermions are intimately connected with the properties of 
the hermitian Wilson Dirac operator Hw for negative mass — M5. Of specific interest is 
the "super-critical" region, —8 < — M5 < 0, in which both Hw and Dw can have zero 
eigenvalues. Th. .eg,™ shows a sur^risiag!^ rich and interesting phase structure, first 



conjectured and described by Aoki [4, 
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58|. 



12l | in both the quenched and dynamical cases 



The Aoki phase was originally defined as a region in which a non-vanishing pionic con- 
densate spontaneously breaks flavor and parity for the two-flavor dynamical case, with an 
associated flavor non-singlet Goldstone pion. In the quenched case, if one examines Green's 
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functions made from a single flavor of quark, (discrete) parity is spontaneously broken. The 
propagator for a flavor- singlet pion describes a massless state at the critical line separat- 
ing the normal from the Aoki phase, but which becomes massive within the interior of the 
Aoki phase. For Green's functions containing two quark flavors, connected by a vanishing 
flavor-breaking mass term, parity and flavor are broken by the pionic condensate ip'y^T°'ijj 
and massless, flavored Goldstone modes exist inside the Aoki phase. 

Since the hermitian Wilson Dirac operator obeys a Banks-Casher relation, a non-zero 
density of near-zero modes is associated with this condensate in both the quenched and 
dynamical cases. This non-zero density is a strong coupling effect, arising from the disorder 
in the gauge flelds which characterizes the strong coupling limit. As the coupling becomes 
weak, the gauge fleld becomes more ordered and such a density of near-zero modes is expected 
to disappear, except for values of M5 where the free Dirac operator itself has zero modes, 
-Mg = 0, -2, -4, -6 and -8. 

However, the picture suggested by the simple summary above is not complete. It was 
observed [12I, [l^ that the interior of the quenched Aoki phase consists of two qualitatively 
different regions bounded by a critical line Qci—M^). Above gc{—M^) near-zero modes exist 
which are delocalized while below gc{—M^), only localized near- zero modes appear. In 
both regions the pionic condensate and density of near- zero modes are non-zero, but their 
localization properties differ fundamentally. 

This picture was substantially refined and solidified by Golterman et al. Q, H, 0, 



17( 1 who, in close analogy with consequences of randomness in condensed matter physics, 
introduced to QCD the concept of a non-zero mobility edge Ac, as the critical eigenvalue 
of Hy/ above which eigenstates are extended and below which all eigenstates are localized. 
They also applied the McKane and Stone localization escape from Goldstone's theorem to 
lattice QCD at non-zero lattice spacing in the quenched case. Specifically, they showed 
that two quenched fiavors of Wilson fermion with a non-zero mobility edge can display a 
spontaneous breaking of a continuous symmetry without a corresponding Goldstone boson. 
This observation is key to the correct functioning of all the various formulations of Kaplan 
fermions wherever the kernel Hy/ has a non-zero density of low modes. 

This qualitative picture is displayed in Figure [H Following conventional terminology, 
we refer to the Aoki phase as that region of this diagram in which {ipj^T^-tp) ^ and 

> g1{—M^). In this region the coupling is sufficiently large that the mobility edge is zero 
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and long-range correlations result from massless flavor-non-singlet pions. The near-zero 
modes of Hw are delocalized. This is a dangerous region for chiral fermions, with lattice ar- 
tifacts producing unphysical, long-distance correlations. For weaker coupling g"^ < g'^{—M^), 
these long-distance correlations have disappeared. The Wilson Dirac eigenvectors with low 
eigenvalues are localized on the lattice scale and all delocalized modes have eigenvalues 
which, if expressed in physical units, are at least as large as Ac(— M5)/a. These delocal- 
ized low modes do not introduce unphysical non-locality into the corresponding formulation 
of lattice chiral fermions but, as discussed in the next section, they enhance the explicit 
violation of chiral symmetry for domain wall fermions with finite L^. 



III. LOW-LYING SPECTRUM OF Hw AND LATTICE CHIRAL FERMIONS 

We now discuss in more detail the relation between the spectrum of the hermitian Wilson 
Dirac operator Hw and the properties of domain wall fermions. Let us first address the issue 
of locality. While the domain wall Dirac operator contains only nearest-neighbor couplings 
in five dimensions (i.e. is ultra- local), the four dimensional effective low energy theory which 
it produces may contain unphysical non-locality. One approach to investigating this question 
is to examine the locality of the overlap operator that results in the Lg ^ oo limit 



18 



19|. 



This was done by Hernandez et al. [20!] using Hw as the kernel in the sign function of 
the overlap operator. They showed that a gap in the spectrum of the eigenstates of Hw 
guarantees locality of the resulting overlap operator. The proof rehes on the ultra-local 
nature of Hw implying ultra- locality of a finite polynomial of Hw- The "Shamir" kernel Ks 
for the overlap operator that corresponds to the = oo limit of domain wall fermions can 



be taken as Ks = Hwi'^ + Dw } 



21 



22 



23| and it is not manifestly ultra-local. However, 



the operator 2 + Dw in the denominator of Ks is the Wilson Dirac operator with a 2 — M5 
mass term and hence, provided M5 < 2 is outside the super-critical region ensuring that 
this kernel is exponentially localized. The Hernandez et al. proof is therefore applicable to 
domain wall fermions under the modification that a finite polynomial of an exponentially 
local operator is also exponentially local. 

However, as discussed above, the assumed gap in the spectrum of Hw is neither expected 
nor observed in lattice QCD, at least for the range of couplings in which calculations are 
normally performed. Fortunately, as argued by Golterman et al. the detailed properties 
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of the Wilson Dirac operator in the super-critical region described above are sufficient to 
imply that the overlap operator is appropriately local (provided we are outside of the Aoki 
phase and Ac(— M5, g) > 0). If, as is the case in this picture, one assumes that the spectrum 
of extended states shows a gap (the region < A < Ac) so that any eigenstates within the 
gap are localized, then the resulting overlap operator will also be localized (even though the 
density of low modes is non-zero). This is supported by numerical evidence at currently 
affordable couplings. 

Finally, we examine a second difficulty faced by domain wall fermions and related ap- 
proaches, which is closely connected with these low-lying modes of Hw- The residual chiral 
symmetry breaking effects seen at finite Lg for domain wall fermions depend in detail on 
the densities and sizes of the modes of the hermitian matrix Hx. This is reviewed in Ap- 
pendix |X1 and Ht is defined in Eq. IA15[ Ht is used to construct the Fock-space transfer 
matrix in the fifth dimension, T defined in Eq. IA13[ This transfer matrix describes explic- 
itly how the left and right walls are coupled when Lg is finite. While the relation between 
the operators Hw and Ht is somewhat complex, it can be shown that their zero modes 



coincide |3|, and also that -logT = Ht = 2 tanh~^ Ks, where Ks = [2l|, [22, l23|. 

Further, the corresponding approximation to the overlap operator, and its — >■ cxo limit 
are: ll 



^approx ^ 1 [1 + ^ + (1 _ tanh(^ifr 
Dov = i [1 + m + (1 - m)75SgnifT)] , 



(1) 



or equivalently, 



(2) 



^approx ^ 1 [1 + ^ + (1 _ tanh(L, tanh"^ 
Dov = |[l + m + (l-m)75Sgni^5)]. 

The equality of the operators sgnHT and sgni^'^ follows easily from the relation Ht = 
2tanh.~^ Ks and the recognition that the inverse hyperbolic tangent is a function which 
preserves the sign of its argument. 

Thus, the near-zero modes of Hyy discussed above are expected to correspond to near- 
zero modes of Ht, modes which (in the absence of an explicit mass term) will dominate the 
coupling of the left- and right-handed sectors of low-energy domain wall QCD. Similarly, the 
mobility edge structure of Hyy is also expected to describe Ht- As is exploited in the next 
section, this implies a simple structure for the asymptotic L^-dependence of the residual 



mass. The near-zero localized modes will give a power behavior (l/L^), while the extended 
modes above the mobility edge give an exponential fall-off {e^^"^" / Ls). This dependence 
will be used as one of several diagnostics to investigate the nature of the Ht spectrum in 
our simulations. 

The strong-coupling behaviour of the Wilson Dirac operator does not leave one at liberty 
to simulate QCD at arbitrarily coarse lattice spacings. For a given lattice action, the phase 
boundary where the gap in the spectrum of extended states vanishes defines the coarsest 
lattice spacing at which the formulation makes sense and implies a minimum cost that must 
be paid to have the formulation under control. The phase boundary is action dependent 
and, in practice, we will wish to stay well within the phase such that the mobility edge is 
not small. 

We have performed the first numerical study of the localization of the hermitian Wilson 
Dirac operator with three mass degenerate flavors of light dynamical quarks and demonstrate 
that a programme of 2+1 flavor dynamical DWF simulations of QCD is rendered affordable 
by our current QCDOC computer systems. These considerations are directly relevant for 
all chirally symmetric lattice QCD formulations, and we shall later illustrate the smooth 
connection between finite Lg and the Lg —>■ oo overlap limit with numerical data for the 



plaquette. Our work has been presented at Lattice 2005 [2J], and later, while this paper 



was being completed, a broadly similar study was presented at Lattice 2006 using dynamical 



overlap fermions 25 1. 



IV. RESIDUAL MASS 

In this section we discuss the residual mass in some detail, including a careful discussion 

urman and Shamir. 



of its dependence on Lg as predicted by the transfer matrix analysis o 
As we will demonstrate and has been worked out previously [,1^, \l5 



16 



17 



24 



26|, this 



dependence on Lg permits the effects of both the localized near-zero modes and the extended 
modes above the mobility edge to be recognized. This gives useful information about the 
general character of chiral symmetry breaking as well as an understanding of the origin of 
the residual mass itself. 

In the limit of small lattice spacing, or equivalently weak gauge coupling, the spectrum 
of the domain wall Dirac operator for a typical gauge background is expected to have a 
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physical, four- dimensional component with eigenvalues A ~ Aqcd as well as unphysical, 
five-dimensional states with A ~ 1/a. The low-lying modes permit the accurate simulation 
of QCD while the large eigenvalues are lattice artifacts similar to, but more numerous than, 
the unphysical large eigenvalues found in other lattice fermion formulations. Ideally, the 
physical four- dimensional states will be bound near the s = and s = — 1 walls. To the 
extent that Lg is large, there should be little mixing between the left-handed chiral states 
localized on the left {s = 0) wall and the right-handed chiral states localized on the right 
{s = Lg — 1) wall. 

The low-energy properties of this approximately chiral theory can be described by an 
effective Lagrangian, C^s- If effects coming from non-zero lattice spacing or arising from the 
overlap between the left- and right-handed states are neglected, this effective Lagrangian 
will be precisely that of QCD. The corrections coming from these effects can be described by 
adding extra operators to Ces- The only relevant operator of this sort is a dimension-3 mass 
term. The next most important term is the familiar dimension-5 Sheikholeslami-Wohlert 
term. Thus, 



In this paper we focus on the coefficient, rrires, of this residual mass term as an important 
measure of the violation of chiral symmetry arising from the finite extent (Lg) in the fifth 
dimension. In this section, we will discuss the dependence of mres on Lg by using the transfer 
matrix formalism of Furman and Shamir j^, described in Appendix 1X1 which is a prerequisite 
for those unfamiliar with this formalism. 

A. Transfer matrix in the fifth dimension 

Using the transfer matrix approach, the approximate chiral symmetry of domain wall 
fermions is realized by writing a general Green's function by semi-independent products of 
left- and right-handed operators. As is described in Appendix |Al if we arrange the field 
variables appearing in a general Green's function by segregating the left- and right-handed 
fields into two products, CL[PLg, gPfj] and OR[PRq,qPL] then, we can write 



Ces = jCqcd + m.esqq + c^qa^T^^q + ... 



(3) 



{OL[PLq,qPR] On[PRq,qPL])L^ 



(4) 




Pj.lf'^^ORlPRd^a^PLPimf)} 
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The left-hand side of this equation is a general, DWF Green's function defined as a Grass- 
mann integral over the 5- dimensional domain wall field variables "^{x, s) and \E'(x, s) using 
the usual DWF action, Sp = — ^Z}^^^\E'. Here the Grassmann fields \E'(x, s) and \E'(a;, s) are 
functions of a space-time coordinate x and fifth-dimensional coordinate s, < s < Ls — I. 
The domain wall Dirac operator D^^^ is that of Shamir and Furman and Shamir [sl 
and in our notation is given by 

D^JJ'A^^^ ^f) = ^s,s'DlAm + (5) 

+ (Ms -4)5,,,. (6) 



2 



(1 - 75)5s+i,s' + (1 + 75)^5-1,5' - 25s,, 

(1 - 75)5s,L,-l5o,s' + (1 + 1^)5s,o5Ls-1,s' ■ (7) 



The physical four-dimensional Grassmann fields, q{x) and are constructed from the 

five-dimensional fields \E'(x, s) and '^{x, s) according to Eqs. IA16l and IA17[ 

The right-hand side of Eq. H] is the trace of a product of operators acting on a many- 
particle Fock space. In contrast to the usual field theory setting, these particles are located 
in a four-dimensional coordinate space. Thus, the fields a, and a\. depend on the space-time 
position X, as well as other spin and color indices which are not shown. These operators, 
as well as the transfer matrix T, and mass operator 0{mf), are defined in Appendix lAl 
The factor Z{Ls) is an Lg-dependent normalization factor. The mass M5 is the domain wall 
height. The mass m/ is the bare mass of the physical fermions. 

In the limit in which mj = and Lg 00, T^" is viewed as projecting onto the state with 
the largest eigenvalue (normalized to be unity), and the right-hand side of Eq. Hlbecomes the 
product of two completely independent matrix elements corresponding to isolated left- and 
right-handed theories (see Eq. IA35I) . This exact chiral symmetry is violated by the operator 
0{mf) which represents the explicit mass term and mixes the left- and right-handed factors. 
Similarly, the corrections to the asymptotic limit T^^ — > |OH)(0/f| of Eq. [A36] introduce chiral 
symmetry breaking arising from finite Lg. 

In the above discussion, and that to follow, we are ignoring the effects of anomalous chiral 
symmetry breaking. Such effects require modifications to our discussion for background 
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gauge fields with non-zero Pontryagin index. For such gauge fields the state \0h) will carrying 
a fermion number different from that of lO^), thereby introducing flavor-singlet chirality 
correlations between the otherwise independent left and right-hand factors in the flrst term 
on the right-hand side of Eq. H] 28|]. In the interests of simplicity, we do not consider such 
gauge conflgurations here. However, the effects of these conflgurations do not alter our 
conclusions. 



B. Residual mass and the transfer matrix 

We will now explicitly study the flrst corrections to the large Lg limit of in Eq. H] 
where, for clarity, we set m/ = and use the relation 0{mf = 0) = |05)(05| discussed in 
Appendix [Aj 

{OL[PLq,qPR]OR[PRq,qPL])j^^ = Z{L,) {Os\OL[PLa, -d^ Pr] (8) 

N+ N- 
k+=l k~=l 

OR[PRa,d^PL]\Qs). 

Here the leading term, the projection operator \()h){^h\-, divides the Green's function into 
two independent factors demonstrating the separate flavored chiral symmetry of the left- and 
right-handed fermions. The flrst correction permits quark number ±1 exchanges between 
these two otherwise independent sectors. This suggests that the second and third terms in 
Eq. [8] should be interpreted at low energies as the residual mass operator nij-esgq expressed 
in this 5-dimensional language. 

A relation between this residual mass term and the second and third terms in the operator 
in curly brackets in Eq. [8] can be obtained if we express the operators 6^ and pk in terms of 
the conventional 4-dimensional fleld d^, inverting Eq. IA25I 

- y^^mresq{Xn)q{Xn) -» ^res | -Pfi)Q I Og) (Og I (-P-RQx„ )a (9) 

n n 

+ {PLd.X\^H){^H\{dl^PL)c] 

^ {dl^PR)M{^H\{PRdyJp 
k n,n' 

k n,n' 
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(10) 

Here the RHS of Eq. [9] is the usual residual mass term, —mj.(,sqq, written in the transfer 
matrix language (see Eq. IA19I) , while the expression in Eq. [10] is a version of the leading 
order contribution from the second and third terms in Eq. [HI The projection operators Pr 
and Pl can be introduced into the right-hand side of Eq. [10] because the other chirality will 
vanish in the limit mj = [59]. 

We can justify this relation and estimate m^es if "we assume that the sums 
Sfc 'i^fe (^n) Vfc (l/n') are localized on the long distance scale at which mres is defined: 

KaA^, y) = J2 'l^tai^y'PkAyy'''^^ ^ rn,,s5\x - y)5^,p , (11) 

k 

where a and (3 are each spin and color indices, and we assume that the diagonal spin/color 
structure will appear and any dependence on the label ± will disappear when a volume 
and/or gauge average is performed. We can then approximate: m^es = / d'^^Ka^x^y) ~ 
Y^R^Ka^yyy), where the "radius" R estimates the small region in x that contributes and 
summation is performed over the repeated spin and color index a. Finally, K{y, y) and 
hence mres can be determined by integrating over y and using the orthonormality of the 
eigenfunctions ^^^(y): 

R^ f 



' k 



poo 

Ve-^^^^* = i?M rfAp(A)e-^^^ (13) 
I, Jo 



RtPeiXcf-^ + RfpiiO)^. (14) 



Here p(A) on the RHS of Eq. [13] is the density of eigenvalues of Ht per unit space-time 
volume, color and spin. The final Eq. [TJ] is a generalization of Eq. [13] displaying the ex- 
pected contributions of extended (e) and localized (/) modes with possibly different radius 
parameters, Ri and Re- Eq. [H] should be viewed as a theoretically motivated model for 
the large Lg behavior of the residual mass where Ri is the radius of the localized states 
with eigenvalue close to zero. The radius associated with the extended states, Rg, arises 
from the limited spatial resolution provided by a sum over extended states resulting from 
the high-momentum cutoff intrinsic to the lattice formulation. This suggests Re ~ a. In 
this final equation we have also taken the limit of large Lg to simplify the integral over A. 
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This is the "standard" result for the dependence of m^es on Lg, although the factor of l/Lg 
is usually missing from the left-most term. We will use Eq. [H] elsewhere in this paper to 
interprete the dependence of m^es on Lg. 



C. Residual mass and the 5-dimensional axial current 

While the above discussion provides a direct connection between corrections to the Lg 
oo limit and the residual mass, it does not give a simple way to compute rrires- This is 
typically done using a partially conserved 5-dimensional axial current introduced by Furman 
and Shamir 

Ki^) = E - (15) 



s=0 

where the current s) is defined by 

jI{x, s) = ^ ^{x + fi, s){l + -ff.)Ul^^j'^{x, s) - ^(x, s){l - -f^p^j'^ix + /i, s)j , (16) 

and 6 is a flavor index and t'' is a generator of the flavor symmetry. The equations of motion 
imply the conservation law 

A.Alix) = 2mfq{x)t'^'q{x) + JI^{x) (17) 

where Jg^ is defined in Eq. IA20I and we define {A^f){x) = f{x) — f{x — fi). 

In the limit of large Lg and small lattice spacing, the extra Jg^ term in the conservation 
law in Eq. [T7] should describe the chiral asymmetry of the leading dimension-3 term in the 
domain wall fermion effective Lagrangian, niresW- 

4^ ^ 2mresgtS'g. (18) 

This equation can be used directly to determine m^es by evaluating the ratio of the matrix 
elements of Jg^ and qt'^'j^q. In this paper, we use the ratio of pion correlation functions: 

{2qt^^^{x) qt^lHiy)) ^ ' 

evaluated for \x — y\ » a to determine m-ees numerically. 

As a consistency check, we should demonstrate that the definitions of mj-es given in Eqs. [10] 
and [18] agree. We begin by generalizing Eq. IA22I to express a general Green's function of 
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J^q{z) with other physical fields, in a form in which the latter are separated into the left- 
and right-handed factors Ol[Pl(1,(iPr] and (^^[PRg, gPi]: 

{Ol [PlQ, qPr] JUz) On [Pnq, qPl] ) 

= Z(L,)(05|OL[PLa,-atp«]f^att''a,f^O«[P«a,atp«]|05). (20) 

Since we wish to compare with Eqs. [9]and[T0|, we can simplify this expression by integrating 
the position variable z over space-time: 

N+ N- 

Z(L,)(05|OL[PLa, -a^PR]f'f[Y. - J] p^.t^p^-jf '^^^[P^a, atp^]|05). (21) 

fc+=i fc-=i 

Here we have also replaced the operators and dj. with those creating eigenstates of T, 
using the definition in Eq. IA25I and the orthogonality relations obeyed by the coefficients 
Next, we can evaluate the factors in the limit of large Lg using Eq. IA36[ The 
leading contribution to Eq. [2T]results when the next leading contribution to (the second 
and third terms in Eq. IA36P is substituted for each factor: 



fc+=l k~=\ 

N+ N' 

{Y. ^l^'^^- - pWp>^-} (22) 

k+=l k-=l 

(5^e" ^t^6lM{0H\d,^+J2^~ ^^^pl.\0H){0H\Pk-)OR[PRa,a^PL]\0s). 



Exploiting the simple form of the operator in the curly brackets, Eq. [22] can be considerably 
simplified, yielding: 

Z{Ls){0s\OL[PLa,-d^PR] (23) 

N+ N- 

1^ e"^^+^=6i+|0H)t'(0i^|6fc+ - Y e-''^-'''pl-\0H)t'{0H\Pk^]OR[PRa,a^PL]\0s) 

k+=l k-=l 

where the fiavor generator acts on the suppressed fiavor indices carried by the operators 
6^+, 6k+, p\- and pk- ■ Now we can directly compare the operator appearing in curly brackets 
in Eq. [22] (which should reduce to the operator 2mres^t^7^q' at low energies) with the second 
and third terms of the operator in curly brackets in Eq. [8] (which we have shown does reduce 
to the operator mres^? in this limit). These two operators are directly related at first order in 
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S by an infinitesmal chiral transformation where operators in the left-hand factor transform 
as 

ol^ ^ + iSt") pl pl_il - i5t') (24) 
while operators in the right-hand factor transform as 

6fc+ ^ (1 + i5t'')dk+ pk- ^ (1 - tSt')Pk- . (25) 

Since this relation between the operators in curly brackets in Eqs. |23] and [8] will continue 
to hold in the long distance limit, each limiting quantity, 2mj.esqt'''y^q and niresQQ, must 
contain the same value of mres- Thus, the residual mass defined in Eq. [19] agrees with that 
appearing in the residual mass operator deduced from the large Lg limit of the transfer 
matrix for the case rrif = 0. This demonstrates the expected consistency with a Symanzik 
effective Lagrangian description of the low-energy effects of mixing between the s = and 
s = Lg — 1 walls. 



V. LATTICE ACTIONS AND ENSEMBLES 



Our simulations have been performed with three flavors of dynamical DWF for a class of 



improved gauge actions. Our notation is the same as in 
here. The partition function we simulate is 
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3l| , which we briefly review 



//* 3 /* 3 

i=i i=i 



(26) 



where the index i runs over the up, down and strange quarks. The Pauli-Villars fields $pv,i 
are needed to cancel the bulk infinity that would be produced by DWF as ^ oo. In 
particular, the total action is 

S = SciU) + Sf{^, U) + Spv{^\ U). (27) 

The gauge actions we consider are of the form 



Sg[U] 



3 



(28) 



;i-8ci) J2 P[Ul,,^u + Ci J2 ^P^^ 

where P[U]x,i_iu and R[U]x,iiu represent the real part of the trace of the path ordered product 
of link variables around the 1x1 plaquette and the 1x2 rectangle, respectively, in the 
/i, z/ plane at the point x, and P = with g the bare quark-gluon coupling. We use two 
common choices for Ci. 
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the Iwasaki action which sets c\ = —0.331 32, l33|]; and 



the DBW2 action which has Ci = -1.4069 M,m\. 



Both are known to reduce the residual mass very effectively in the quenched case. 

The fermion action in Eq. [27] can be expressed in a form very close to the way it appears 
in the 3-fiavor numerical work if it is written as: 

3 



1=1 



1/2 



(29) 



where is the input bare quark mass for the ith light quark flavor and the DWF Dirac 
operator -Ddwf is deflned in Eq. |5l We only consider the case where all light quarks have the 
same value for the flve-dimensional domain- wall height, M5. The action for the Pauli-Villars 
fields is similar, except that the quark mass is replaced by 1, to yield 

3 



Spy = pi3WF(^5, 1)/^DWf(M5, 1) 



1/2 



(30) 



i=l 



The resulting determinants are suitable for simulation using a weight 



n 



det^/2 






Dj3WF(^5,l)^DWF(M5,l) 





(31) 



The action may be simulated using either the inexact R algorithm 36l|, with associated 



step size errors, or the exact RHMC algorithm of Clark and Kennedy |37l . I38l ]. For the 
RHMC simulations, each of the square root factors in Eq. [3T]are approximated by a rational 
function accurate over some eigenvalue range. Thus our simulation uses a 1+1+1 framework 
that happens to be mass degenerate. Where our mapping of the Aoki phase finds a clearly 
non-zero mobility edge it should be relevant to our more general programme of 2+1 fiavour 
simulations where a chiral limit is taken in mass degenerate up/down quark masses at 
fixed strange mass. When written as a partial fraction expansion, each term is individually 
amenable to the usual HMC pseudo-fermionic integral (though, of course, a single pseudo- 
fermion and multi-mass inverter is used for all terms in the partial fraction expansion). We 
used the R algorithm for the earlier runs in our parameter search and used the RHMC out 
of preference, as documented in Table [B Our implementation of the RHMC used separate 
fields as stochastic estimators of the numerator and denominator of Eq. [3T] We simulated 



all three fiavors in the numerator as separate single-fiavor species, and treated the Pauli 
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Villar fields as separately estimated two- and single-flavour pseudofermions. Thus, we did 
not make any explicit use of the mass degeneracy between up and down quarks, (as would 
have been the case, for example, in two-flavor HMC). 
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m 



4l| , and have 



The ensembles used in this study are a subset of those presented in 
been used for precision valence measurements with Lg = 8. We used a reduced statistic 
subset to perform a consistent study as we vary Lg in the (still ex pen sive) valence analysis. 



The larger statistical sample used for the Lg = 8 data points in [39|, |40|, |4l| lead to lower 
errors on results for mj-es than we present in this work. We include only lower precision 
values in this paper to keep the = 8 analysis identical in all systematic respects to the 
Lg > 8 points that are unique to this work. The ensembles used in this work are listed in 



39, 



m 



4l| where the entire 



Table [H For the full simulation details, we refer the reader to 
dataset is presented. 

All ensemble and correlation function computations were performed using the Columbia 
Physics System (CPS) on QCDOC computers in Edinburgh and Brookhaven. All measure- 
ments of four-dimensional eigenmodes were performed using the CHROMA physics system 
42], while the 5-dimensional modes were measured using CPS. The relation between these 
modes provides a non-trivial cross-check between code bases. Both code bases made use 
of the optimized Bagel QCD library {4^. Visualization was carried out using the OpenDX 
software package. 



VI. SPECTRUM OF THE HERMITIAN WILSON DIRAC OPERATOR 

In this section we study the eigenspectrum and eigenmodes of Hw on our ensembles. We 
determined up to 256 eigenvectors of H\y on 25 configurations per ensemble, at the negative 
mass M5 = 1.8 that was used in our domain wall action for the sea quarks. Additional 
eigenvalues were generated on a number of configurations as a function of M5. Our study 
therefore includes the dependence of the low- lying spectrum on M5, often called spectral 
fiows, in Section IVIAi We also study the spectral density of near-zero modes in Section IVIBl 
the microscopic shape distribution of the low-lying modes in Section fVI CI and the relation 
between modes of the hermitian 5-dimensional operator i/owF = Ts-Rs-Ddwf and Hw in 
Section IVID[ where is the refiection operator in the fifth dimension. 
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A. Spectral flows 



Prior to recent work on localization, it was believed that a gap in the spectrum was 
required for the correct behaviour of the various lattice formulations of chiral symmetry. In 
this context, it was believed that by obtaining spectral flows of the lowest few modes on 
a "typical" configuration, the health of the formulation was indicated by a nicely opened 
"eye" for some range of M5 containing the Wilson mass used in the simulation. Recent 
improved understanding suggests that a non-zero, but low, density of localized states exists 
throughout the region in which we can afford to simulate. A gap, as such, does not exist 
after an ensemble average and, if spectral flows from sufficiently many configurations were 
overlaid, the entire eye would be filled in. 

However, on a small sample of configurations an open eye is indicative of a low density 
of low-lying modes and, we argue, a non-zero mobility edge since our results show that 
the spectral density grows rapidly across this edge. The flows presented here also allow 
comparison between our ensembles and those of previous works. These plots for the lowest 
modes. Figures [2] and [31 are therefore useful information at the larger values of (3 where 
the open eye is indicative of a non-zero lower bound for the mobility edge, as presented 
previously [44 1. 

Our spectral flow results are certainly indicative of a non-zero mobility edge for (3 > 0.78 
for the DBW2 gauge action and (3 > 2.2 for the Iwasaki gauge action. For smaller values of 
f3 the spectral flows are inconclusive as no opening is displayed suggesting a non-zero density 
of low modes. We shall return to look at the localization structure in a more sophisticated 
fashion and classify these modes in Section IVI C[ 



B. Spectral density 

For each dataset comprising nconf configurations we label the ith eigenvalue on the jth 
configuration Xij. We use a somewhat simplistic binning approach to estimate the spectral 
density for a partitioning of the spectrum. For the kth bin we estimate the density as 

P ([-^min + ^max]/2) = ~ TT _ ) 

"conf ^max ^min 

where Ni^ is the number of eigenvalues X^j G [A^m + -^max] 5 i = 1, . . . ,nconf- Figure H] displays 
our results for the spectral densities for our various datasets, using a fixed bin width of 0.02. 
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Some clear, and unsurprising, trends are evident: the spectral density p(0) is non-zero, and 
falls rapidly as (3 is increased. The quantity p(0) is slightly higher for the Iwasaki action at 
(5 = 2.2 than for DBW2 at /? = 0.78 which have comparable lattice spacings. The density 
p(A) rises rapidly, and perhaps exponentially, in A as we move away from zero. 



C. Microscopic study of the mobility edge 

We wish to study the size distribution of our eigenvectors, ipi^x), as a function of the 
eigenvalue in order to test the mobility edge conjecture. This approach is novel and was 
first presented in 2J] for dynamical DWF, with a similar approach taken in 25| for dynamical 



overlap. A standard measure is the inverse participation ratio. 



A normalized extended state will produce P ~ V^, while a completely localized state will 
produce P 1. Thus the (inverse) participation ratio essentially counts the number of 
occupied sites and, in the infinite volume limit, is a useful order parameter for delocalization 
transitions. However, this metric is not sensitive to the shape of a state and so makes no 
distinction between occupying two adjacent sites, and occupying two well separated sites, 
such as the illustrative non-local example in Figure [5l 

We want to understand how to reliably distinguish exponentially compact states from 
those with any degree of long-range support. In the infinite volume limit the use of the in- 
verse participation ratio only identifies long-range correlations arising from structures whose 
four dimensional volume diverges and we seek a more robust measure. 

While a better measure could be a moment of the eigenvector density, such as 



x)\x\ 



we have found empirically, by visualizing data, that near the delocalization threshold the 
dominant long range correlation arises from multi-peaked eigenvectors, which fall exponen- 
tially around multiple centers. We define a more robust measure of exponential localization 
length, which generalizes the one dimensional case of Figure |5]to higher dimensions as follows 

• The mode density is defined as = ip^ {x)4'{x). 
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• We identify a center xq such that 0(xo) > (j){x)Wx ^ xq. 

• For every site, we define a radius r{x) = \x — Xq\ using the periodic mirror image 
nearest to xq. 

• We define an effective locahzation exponent L^six) = iog(i>{xo)~iog(j>{x) ■ 

• Finally, we define a robust localization length L^ax = max Leff(x). 



As the eigenvalue approaches the mobility edge from below, we observed two clear pro- 
cesses by which delocalization takes place. Firstly, the mean rate of fall off from the center 
decreases creating much bigger states. Secondly, states become increasingly multi-centered, 
with exponential fall off between a number of satellite peaks. The above definition of lo- 
calization length is designed to refiect both the process of single-peak broadening, and the 
development of multiple peaks. The inverse participation ratio would be a less robust mea- 
sure of locality in that it would not differentiate two well-separated 6 functions from those 
on two neighbouring sites, or indeed would not identify a lower-dimensional extended sheet 
in the infinite volume. We display scatter plots of this definition of localization length for 
the 256 lowest eigenvectors of on 25 configurations per ensemble in Figure [6] and [7] for 
the DBW2 and Iwasaki gauge actions at /? = 2.13 and f3 = 0.764 respectively. 

D. Relation between eigenmodes of the five and four-dimensional operators 

We computed the lowest 10 modes of the 5-dimensional hermitian domain wall operator 
Hbwf = 75-R5-DDWF, using the Columbia Physics System, on a single configuration from a 
16^ X 32 X 8 DBW2 three fiavor ensemble at /3 = 0.764. On the same configuration, we used 
the lowest 256 eigenmodes ipj{x) of Hw, computed using the CHROMA software package, 
to express the 5-dimensional modes s) as a sum over 4-dimensional modes using 



r(z)>5 




(32) 



where 



N^ = J2'^lix,s)^i{x,s) 



X 
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normalizes each s-slice of the 5-dimensional eigenvector to unity. This means that if the 
basis of 4-dimensional eigenvectors were complete 



Y.a.Asf = l. (33) 



However, since we are only using the lowest 256 eigenmodes, we expect that Eq. [32] will, 
at best, be true only deep in the bulk for large Lg where low modes of the transfer matrix 
dominate, and only up to differences between Ht and H\y. 

We therefore computed 5-dimensional modes for a valence Ls = 16 (N.B. not the unitary 
Ls = 8 case) finding that, while the basis is very much incomplete on the wall, the description 
it gives in the 5-dimensional bulk is as much as 80% complete by the mid-point, despite the 
mismatch between Ht and Hw- (This also serves as a useful cross check between the 
eigensolvers in the two independent code bases.) 

The coefficients \aij\'^ for a typical 5-dimensional chiral mode are displayed in Figure [TOl 
and it can be seen that only 3 low-lying modes of Hw dominate the s-dependence of eigen- 
modes of i^DWF- Thus, we clearly demonstrate from numerical data that a few localized low- 
lying 4-dimensional modes of H^ almost completely describe the coupling between walls for 
each 5-dimensional mode of the hermitian DWF operator for large Lg. As discussed, these 
modes occur with a low density - perhaps a "dilute gas" is an appropriate picture - but, for 
sufficiently large Lg, they dominate the contributions to chiral symmetry breaking due to 
exponential suppression of the extended modes above the mobility edge, see Eq. [HI 

The localized low modes of Hw thus play a role in providing conduits into the bulk for 
chiral symmetry breaking, causing corresponding localized spikes in the correlation function 
for TTires- This is corroborated by the iso-surfaces (produced with the OpenDX package) 
of constant ifowF-eigenvector density (i.e. surfaces defined by \I'^(x, s)\l/(x, s) = constant) 
obtained in our simulation. Examples of this are shown in Figures [8] and [9], which represent 
lower dimensional views, at different spatial locations, of the same eigenvector. A movie, 



60|. 



raster scanning x and ?/ as a unified "movie-time" , can be obtained from 

Since the low modes of Hw are scarce and very localized, it is to be expected that if we 
are in a region of parameter space where their contribution is the dominant contribution to 
chiral symmetry breaking, quantities such as mres will display larger non-Gaussian tails (or 
"lumpiness") than if the exponentially falling contribution from extended states is dominant. 
Thus, one should not be surprised by the need for very high statistics to achieve a convincing 
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estimate for mres in this region of parameter space. 



VII. RESULTS FOR THE RESIDUAL MASS 

As discussed in Section llVt we calculate the residual mass from a ratio of correlation 
functions containing the axial Ward identity defect, Jg^, defined in Eq. IA20I 

where 6 is a fixed fiavor index, and J^{y,t) is the local 4-dimensional pseudoscalar density 
defined on the s = and s = Lg — l walls. The operator P^(0) is either J^{y, t) or a spatially- 
smeared variant defined at t = 0. The residual mass is determined from a constant fit to 
R{t) over a range of times, t, that are sufficiently large that lattice artifacts should be absent. 
Recall that the ratio R{t) should equal the residual mass provided that the numerator and 
denominator of Eq. [3l] are dominated by physical states. There is not a requirement that 
only a specific lowest state contributes. For the calculations reported here this average is 
performed over the range 7 < t < 16. Figure [TT] through Figure [TSl contain plots of the fits of 
R{t) to a constant. We use a combination of point. Coulomb gauge-fixed Wall and Coulomb 
gauge-fixed exponential sources on the various datasets listed in Table [Til and display the 
fitted values for m^-es in Table [ITT] and Table [IV] for the DBW2 and Iwasaki gauge actions 
respectively. 

In order to conserve computer resources, we studied the dependence of mres on Lg by 
working with our fixed set of gauge configurations generated with = 8 but varying 
the "valence" value of Lg that appears in the propagators used to compute the ratio in 
Eq. [Ml The resulting Lg dependence of the resulting m^es can be viewed as a (likely good) 
approximation to what would have been obtained if we also varied the Lg appearing in the 
fermion determinant. This also can be viewed as a self-consistent study the properties of 
the 5-dimensional transfer matrix defined on gauge configurations generated with this fixed, 
Lg = 8 fermion action. 

Next we examine the degree to which Eq. [TH describes the dependence of m-res on Lg in 
our simulations. Here we rewrite that equation in the form: 

m,,,(Lg) ~ (cie-^=^= + C2) ^. (35) 
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Recall that the first term is expected to come from extended states with eigenvalues A near 
the mobility edge: Xc < X < Xc + j-- This term is exponentially suppressed for large Lg 
and thus can be easily reduced by increasing Lg. The second term arises from low-lying 
localized states with < A < y-. Here C2 is proportional to the density of near-zero modes, 
p(0). We require a demonstrably non-zero mobility edge Ac for safe QCD simulations with 
chiral formulations. It is acceptable to have a significant non-exponential component in 
rrij-es- However, the presence of such a term makes it difficult to significantly suppress chiral 
symmetry breaking effects by simply increasing Lg. 

We show fits to the functional form in Eq. [35] for degenerate three fiavor ensembles at 
several gauge couplings for the Iwasaki and DBW2 gauge actions in Figures [11] and [2D1 The 
fitted parameters are listed in Table [V] While the exponent Ac increases only slowly, the 
parameter C2 falls by orders of magnitude as the lattice spacing decreases. Interpolating 
between the values of lattice spacing given in these tables suggests the localized near-zero 
modes are suppressed more by the DBW2 than the Iwasaki gauge action. This behavior will 
be correlated with the rate of topology change discussed in the next section. 

We note that while the residual mass term may be incorporated in a renormalized quark 
mass, one should extrapolate to = cxd to remove the residual chiral symmetry breaking 
effects of the C5 term in Eq. [3] However, even for finite Lg, such terms are expected to be 
very small compared with even the errors obtainable with non-perturbative improvement 
techniques for standard clover fermions, and thus can be treated as zero. Experience from 
clover fermions suggests that this level of error in the coefficient C5 will certainly leave the 
results unmodified in practice, and that the unsuppressed O(a^) errors will be dominant in a 
continuum extrapolation. Figure [21] displays the dependence of the simplest of observables, 
the plaquette, on Lg for the sea quark content. It can be seen to be very weak (particularly 
for Lg > 16) and consistent with an exponential approach to the ^ 00 limit, and is 
further evidence that all corrections to this limit are exponentially small. 

VIII. LOCALIZED MODES AND TOPOLOGY CHANGE 

Topology change, as defined by the index of an overlap operator will, be accompanied by 
a change of sign of a mode of Hwi—M^) in the region —2 < — M5 < 0. There will be an 
instant r = tqIM^) in the molecular dynamics time at which Hwi—M^) has an exact zero 
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mode for any appropriate choice of M5. Such a level crossing leaves the topological index 
of an overlap operator based on this kernel indefinite, depending on whether — M5 is placed 
to the left or right of the level crossing, and it is reasonable to conjecture that low modes of 
Hw{—M^) play a significant role in topology change in dynamical simulations with current 
algorithms. 

In this section we shall present results for the topological charge time history in our 
ensembles, and consider the relationship between tunneling rates, the eigenmode spectrum, 
and details of the algorithm and action. 

In Figures [22] and [23] we display the topological charge time histories for all our datasets 
measured usingthe same 0(a^)-improved definition of the topological charge density as in 



references 
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461], implemented using the Columbia Physics System. It is certainly the case 



that the tunneling rate is substantially lower in those datasets with lower spectral densities. 
These are also the configurations with the finest lattice spacings. 

Sampling all topological sectors has proven problematic with RG improved gauge actions 
near the continuum limit even in the quenched approximation; some of the authors found 
previously that the DBW2 gauge action samples topological sectors rather poorly at a^^ = 



3GeV in the quenched case 



30|, and this loss of tunneling is likely driven by the gauge 



action. Going nearer to the continuum limit yields an increasing potential barrier between 
topological sectors, arising from the gauge action, which prevents topology change in those 
algorithms that either use a local update, or, even worse, Hamiltonian evolution to propose 
a next configuration. 

Whenever the determinant of Hw does not directly enter the probability weight, the spec- 
tral density of low modes of Hw found for the resulting gauge configurations is principally 
determined by the degree of roughness admitted by the gauge action, likely corresponding to 



fluxon-like configurations 47|, [48|]. We conjecture that for domain wall fermions at practical 



values for Ls the tunneling rate is controlled by this low mode density. 



A. Tunneling with dynamical overlap fermions 

Additional problems have been seen with the introduction of exactly chiral fermions in 
dynamical simulations. Exact implementations of the overlap operator introduce a discon- 
tinuous step in the action when a low eigenvalue of Hw changes sign (and hence when the 
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topological index of the configuration changes). This infinitely narrow step in the action is 
guaranteed to be unresolved by numerical integration for any non-zero timestep, and leads 
to a high pro bability of rejection of the proposed configuration. The reflection/refraction 
algorithm [4^, [sol is one sensible, but expensive, response to this. 



Another response has been to deve 
pseudofermion determinant estimation 
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actions which suppress low modes of Hw with 
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53|, and so suppress topology change. While 



many hadronic quantities are not likely to be sensitive to global topology (especially in the 
limit of large volume), there remains a serious concern that a modification which suppresses 
the change of global topology may also interfere with the creation of a proper distribution 
of local topological fluctuations. In the language of instantons, one should demonstrate that 
such an algorithm leads to a physical density of instanton-anti-instanton pairs - a quantity 
that may be important for QCD even when the global topology is fixed. However, we have 
poor tools for knowing whether this globally non-ergodic algorithm is correctly sampling 
local topological fluctuations, given that we have thrown away our best diagnostic (and 
indeed the indicator that originally flagged the problem), i.e. global topology change. 

So far with DWF we have obtained acceptable tunneling rates, and this is likely because 
the response to the lowest eigenvalues of the tanh approximation to the sign function, Eq. [U 
is spread out at the level of our integration timestep. More specifically the determinant 
whose value is estimated contains a contribution tanhL^ tanh~^ A, in place of sgnA, and the 
transition has a width of 0(;^). If we consider the low mode to be evolving in molecular 
dynamics time, with timestep 5t and at a rate A, then the molecular dynamics integral will 
be accurate provided 

5r « X (36) 

There will be no practical impact on the tunneling rate from increasing Ls while the 
condition Eq. [36] is held, and the problem may, in fact, never arise for any practical DWF 
simulation. In practice, the tunneling rate is determined by other effects, such as the low 
mode density determined by roughness admitted by the gauge action. We have also found 
a significant improvement in the tunneling rate, subsequent to this simulation arising from 
improving the quality of our stochastic estimate of the fermion force [5^, where the use of 
a single stochastic estimate for the ratio of Pauli-Villars and light mass determinants was 
found to allow increased molecular dynamics step size and to increase the tunneling rate. 
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This also suggests a third possible response to this problem that arises in simulations with 
better approximations to the overlap operator than our DWF simulations. That is to relax 
the approximation in the molecular dynamics evolution such that Eq. [2S] is satisfied, while 
keeping the approximation accurate for the accept/reject step. This differential treatment 
of Metropolis and molecular dynamics steps is not possible with DWF simulations due to 
the fully five-dimensional pseudofermion fields, but it is possible for many approaches to the 
overlap operator. 



IX. CONCLUSION 



When constrained to keep Ls affordable, there is some degree of trade-off between the ex- 
actness of chirality and the thoroughness of topological sampling with domain wall fermions. 
Both are determined by the density of modes in the region where the tanh approximation 
is inaccurate, and this density is controlled by either the choice of the gauge action, or the 
inclusion of additional unphysical terms in the action. 

One might pause for a moment to ponder where to compromise. Using computer preci- 
sion, {e.g. single precision) as the standard for achieving chiral symmetry invokes computer 
technology rather than physics for an answer to this important question. The authors' con- 
clusion from this detailed study is that we should require adequate chirality, and by that we 
mean that the average symmetry breaking should not compromise the attractive simplifica- 
tion of electroweak current structure and off-shell 0(a) improvement for renormalization. 

We have demonstrated in this paper that we can obtain rrij-es — 0(10^'^) with both the 
DBW2 and Iwasaki gauge actions in the region of ~ 1.6GeV. We have observed that 
dynamical simulations with the DBW2 gauge action rapidly lose topological tunneling as (3 
is increased. In contrast, headroom is left for going to weaker coupling, while maintaining 
good topological sampling, with the Iwasaki gauge action. This approach is certainly at best 
a stay of execution, as topological tunneling will be lost on sufficiently fine lattices. These 
are inaccessible with our current computers and our conclusion is thus clear - the RBC and 



UKQCD collaborations are running SJ] dynamical domain wall fermions with the Iwasaki 
gauge action and Ls > 16. We started from a^^ ~ 1.6GeV and thus we are able to maintain 
both excellent topological sampling and acceptable chiral symmetry violation. 



26 



Acknowledgements 



We thank Dong Chen, Cahn Cristian, Zhihua Dong, Alan Gara, Andrew Jackson, 
Changhoan Kim, Ludmila Levkova, Xiaodong Liao, Guofeng Liu, Konstantin Petrov and 
Tilo Wettig for developing with us the QCDOC machine and its software. This develop- 
ment and the resulting computer equipment used in this calculation were funded by the U.S. 
DOE grant DE-FG02-92ER40699, PPARC JIF grant PPA/J/S/1998/00756 and by RIKEN. 
This work was supported by DOE grants DE-FG02-92ER40699 and DE-AC02-98CH10886 
and PPARC grants PPA/G/O/2002/00465, PP/D000238/1 and PP/C504386/1. AH is sup- 
ported by the UK Royal Society. We thank BNL, EPCC, RIKEN, and the U.S. DOE for 
supporting the computing facilities essential for the completion of this work. 



APPENDIX A: DOMAIN WALL FERMION TRANSFER MATRIX 



Here we collect some useful formulae that describe the transfer matrix formalism as 
applied to domain wall fermions by Furman and Shamir While we do not present a 
derivation of these formulae, they follow reasonably directly from the results presented by 
Furman and Shamir, Narayanan and Neuberger 55] and Liischer 56 1. 

We begin with the five-dimensional domain wall fermion Dirac operator defined in Eqs. 
[6] and [3 (This definition of the domain wall fermion operator is consistent with our earlier 
papers, but differs by a hermitian conjugate from the original notation of Furman and 
Shamir.) We also identify two important 4-dimensional operators which play a central role 
in Liischer's original transfer matrix construction 



(5 - M,)6^ 
1 



(Al) 
(A2) 



Here we have adopted a spinor basis in which 7^ is diagonal. In this basis the 7 matrices 
can be written 



r 



1 a''] 




fi \ 






10 -il 







(A3) 



and the four 2x2 matrices a'^ appearing in both Eqs. IA2I and IA3I are given by cr'^ = {i, a), 
where a represents the three standard Pauli matrices and |(1±7^) projects onto right-/left- 
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handed states. Thus, the matrices B and C have a 2 x 2- dimensional spin structure and 
appear directly in the matrix dI^,{M5) of Eq. [61 




D'W'h) = I - \^ \ ■ (A4) 

Next we need to introduce three pairs of related fermion fields. The first two pairs are 
the 5-dimensional Grassmann variables '$x,s and '$x,s and the corresponding Fock-space 
operators ip^ and The field operators ipx and can be written in terms of two- 
component field operators and dx and their hermitian conjugates: 

= ct5-5 ) . (A6) 

Here the operators c and d both act on vectors in the Fock space and carry space-time, color 
and spinor indices. As Fock space operators, they can be thought of as 2^^^^^' dimensional 
square matrices if we are working with a. x Lt space-time lattice. We treat their space- 
time, color and spin indices as the index on a QL^Lt x 1 column vector. It is on this vector 
that the QL^Lt x QL^Lt matrix acts in Eqs. IA5l and lA6[ We use the conventions that the 
hermitian conjugate, represented by the superscript f takes the hermitian conjugate of both 
the Fock-space operator and the space-time-fiavor-spin matrix. The transpose, represented 
by the superscript t, acts only on the space-time-fiavor-spin matrix. 
Finally we introduce two further operators: 

R{c,d) = exp|c/*(S-5CS-5)c| (A7) 
W{c\d\cJ) = exp|-cMn(5)c-c/^ln(5)rf}. (A8) 

Using these operators we can then connect the Feynman path integral, which provides 
the original definition of domain wall fermion lattice theory, with a Fock-space operator 
expression. For simplicity, we consider an example of the 5-dimensional fermion propagator, 
because a general fermionic Green's function follows exactly the same pattern: 

/ l[d[^x',s']d[^.',s']e^''°'^^'' 
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'^x„s,^x,,s, (A9) 



where 



w^Ri^^^R^w^iw^M^w^y^ ^C*!^/)! (AlO) 

^ilTT'f) = n|(Ca4 + mfclca)idadi + TUfdlda)^ (All) 



and a represents space-time, spin and color indices. While a complete derivation of the 
result in Eq. lAlOl is beyond the scope of this paper, one of the intermediate steps is shown 
in Eq. IA9I In that equation the Grassmann fields ^ and \l/ have been expressed in terms of 
their chiral components according to: 



VP = I - I M/ = ( # ^« ) . (A12) 



A complete derivation of Eq. lAlOl can be constructed from Refs. [s], [ss], l56| . Note, the Grass 



mann fields and appear in the integrand in Eq. IA9l but are not among the variables 
of integration. Instead, these fields are to be evaluated as —nif times the corresponding 
fields at s = 0: — mj\l/^^Q and — m^vI'f'^Q, respectively. 

Following Furman and Shamir, we identify the Fock-space transfer matrix T as 

T = {det B)WmR^W^ =e-'^^^^^ (A13) 
(d^Y 



where d 



c 



i^=(^d\ ct) and (A14) 



, B + CB-^C^ CB-^ , 
e-^^ = . A15 

The unusual arrangement of the operators ipxi and ip^^ and two missing factors of T in 
Eq. lAlOl can be understood if, following Furman and Shamir, we recognize that the upper 
component of i^ and the lower component of ip commute with R\ while the lower component 
of ip and the upper component of ip commute with R. This allows the operators ip and ip to 
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be extracted from between the factors W^R and R^W^, so that the two missing factors of 
T can be assembled. (Note, the operators ip and do not commute with the operator W^, 
but instead have the non-local factor of appearing in Eqs. IA5I and IA6I removed by the 
proposed interchange). 

We will demonstrate this re- arrangement of operators by working out two examples. In 
the first, we translate the expectation value {q{x)q{x)) (the chiral condensate) into this 
transfer matrix language. Recall that the physical, 4-dimensional Grassmann fields are 
defined by: 

qix) = Pl^,,o + Pi?^.,L.-i (A16) 
q{x) = ^.,oPr + ^.,l.-iPl, (A17) 

where we have used Pr/l = (1 ± 7^)/2- If these expressions are substituted into Eq. lAlOl 
we obtain: 

{qix)q{x)) = -{det B)HT!^W^R{PLi^^)c,R^W^ (A18) 
T^'-^W^RitxPL)aR^W^O{mf) 
+ {detB)^ni!^W^R{1pPR)aR^W^ 

T^^-^W'^R{PRilj^)aR^W^O{mf) 



The final step involves commuting the Pli/j and i/jPr past the factor W^R to stand on the 
far left. Similarly, the operators PriP and i^Pl on the right side can be moved past the 
factor R^W^ to stand directly to the left of the operator 0{mf). The resulting expression 
can be written simply in terms of the operators a and aP: 

{q{x)q{x)) = -(det5)^nr|[(Pia,.)aT^l4^L)a (A19) 

+ {alPR)^T''^iPRa,)a\o{mf) 



As a second example, consider the Green's function {Piq^x) J^g{z) q{y)PL) containing the 
midpoint operator J^^ that occurs in the divergence of the 5-dimensional flavor non-singlet 
axial current introduced by Furman and Shamir [^: 

J',^{x) = -^{x, Ls/2 - l)PLt'^{x, Ls/2) + ^{x, Ls/2)PRt'^{x, L,/2 - 1). (A20) 
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Here 6 is a flavor index and is a flavor generator. 



Using Eq. lAlOl for tliis product of four Grassmann fields we find: 

{PLq{x) Jl^{z)q{y)PL) = -(det 5)^nr|#^i?(Pi4)i?^W^^T^=/'"' 

T^^''-^w'^R(^yPR)R^w'^0{mf)\ . 



(A21) 



Performing the same operator manipulations as used to obtain Eq. IA18I we find: 



{PLq{x) Jl^{z)q{y)PL) = (det 5)^=tr|p,.a,T'^ att^a^T^aJP^OK) } 



(A22) 



Here the unit displacement in s between the two factors in J^^ has disappeared leaving a 
very symmetrical expression. The missing 7^ matrix, which one would expect in an axial 
matrix element, has its origins in the minus sign that is present in Eq. IA6I but absent in 
Eq. IA14[ This factor could be restored if we introduced = a^7^, a 5-dimensional analogue 
of the role played by 7° in the standard treatment of the 4-dimensional Dirac operator. 

As a final topic for this appendix, we will examine the effects of the matrix T^" in Eq. IA18[ 
Again following Furman and Shamir, we introduce the eigenvectors and eigenvalues of the 
operator Ht defined in Eq. IA15I 



where E"^ > 0, and write the operator in terms of this basis, filling the Dirac sea in the 
usual way : 



Let \0h) be the Fock space state \0h) which is the eigenstate of T with the largest eigenvalue. 
This state obeys: 



Ht^.^ = -E,.r,- i<k-<N- 



(A23) 



(A24) 




(A25) 



(A26) 



for all and k and 




(A27) 
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Finally we define the normalized operator T = T/Amax making unity the largest eigenvalue 
of T and \0h) the corresponding eigenstate. The matrix T can be written: 

f = T/A^ax = exp{-a^Ha} /A^^x (A28) 

{N+ N- \ 

- KA+''k+ + Yl ^k-Pk-Pl- \ /^max (A29) 

k+=l k-=l ) 

{N+ N- ^ 

- E KA^o,, - Y E^-Pl-Pk' ■ (A30) 
fc+=l k-=l ) 

Thus, in the limit of large L^, T^" becomes the projection operator onto the state |0h), 
separating the right and left operators in an equation such as Eq. IA18I into two separate 
factors, linked by the operator 0{mf). If we further take the limit mj = and define a 
surface vacuum state lO^), which is annihilated by the operators Cx and dx, 

Cx\Os) = dx\Os) = 0, (A31) 

then the operator 0{mf = 0), given in Eq. lAlll reduces to the projection operator |05')(05'| 
and the left- and right-handed sectors become completely independent matrix elements 
between \0s) and \0h)- 

This can be summarized by evaluating a general Green's function depending on the 4- 
dimensional fields q and q. This is most easily done if the fermion fields in the original 
Green's function are ordered with the left-handed chiral fields {Piq and qPn) on the left and 
the right-handed chiral fields {Puq and gPi) on the right: 

(deti? Vax)^- ^^""^ ^^^^^ 

= jy^l: tT{OL[PLd,-d^PR]T'^^OR[Pnd,d^PL]0{mf)} (A33) 

tT{OL[PLd,-d^PR]\0H){0H\OR[PRd,d^PL]O{mj)} (A34) 

^ {0s\OL[PLd,-d^PR]\0H){0H\OR[PRd,d^PL]\0s), (A35) 

where the quantities Ol and Or are two ordered polynomials in their two arguments. 

Corrections to the limit Ls oo can be written by including eigenvectors of T with 
smaller eigenvalues. Such corrections can be easily obtained from Eq. IA30[ 

N+ N~ 

~ |0h)(0//|+ 5^e-^.'"+^»6[+|0H)(0^^|6fc++ e-^.'-^^p^, |0^^) (0^,1^^- + . . . . (A36) 

k+=l k-=l 
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TABLE I: Ensembles used in this study. We also document the trajectories used for mj-es measure- 
ments, using every tenth trajectory in the given range. The trajectory length was r = 0.5 units of 
MD time. All ensembles used three degenerate DWF flavors of mass m = 0.04, and M5 = 1.8. 



Gauge Action 


P 


Algorithm 


Trajectories 


DBW2 


0.72 


R 


800-1500 


DBW2 


0.764 


R 


800-1500 


DBW2 


0.78 


R 


1000-1600 


DBW2 


0.80 


R 


900-1800 


DBW2 


0.88 


R 


1000-3000 


Iwasaki 


2.13 


RHMC 


1500-2400 


Iwasaki 


2.2 


RHMC 


800-1500 


Iwasaki 


2.3 


RHMC 


800-1500 
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TABLE II: Correlation functions generated for valence quarks with different extents of the fifth 
dimension on each of our Lg = 8 dynamical ensembles. The valence and sea quark masses were 
mud = rug = 0.04. We denote point source correlators as "P" , wall source correlators as "W" , and 
exponentially smeared correlators as "S". 



Dataset 


Ls = 8 


Ls = 12 Ls = 16 Lg = 20 Ls = 24 Lg 


= 


28 Lg = 32 


Iwasaki ^ = 2.13 


S 


S 


s 


s 


s 


- 


S 


Iwasaki (3 = 2.2 


W,P 


P 


p 


p 


p 




P 


Iwasaki /? = 2.3 


W,P 


P 


p 


p 


p 




P 


DBW2 /? = 0.72 


w 


W 


w 




w 




w 


DBW2 13 = 0.764 


w 


P 


p 


p 


p 


p 


p 


DBW2 /? = 0.78 


w 


W,P 


W,P 


W,P 


W,P 


p 


p 


DBW2 /? = 0.80 


w 


p 


p 


p 


p 


p 


p 


DBW2 /? = 0.88 


p 


W,P 


W,P 


W,P 


W,P 




W,P 
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TABLE III: Fitted residual mass values for valence quarks with different extents of the fifth dimension for the DBW2 gauge action on each 

of our Lg = 8 dynamical ensembles. The valence and sea quark masses were m„rf = = 0.04. 



Dataset 


Ls = 8 


Ls = 12 


Ls = 16 


Ls = 20 


i^. = 24 


Ls = 28 


Ls = 32 


/? = 0.72 


1.13(1) X 10-2 


4.4(1) X 10-3 


2.2(1) X 10-3 




1.04(5) X 10-3 




7.0(5) X 10-"^ 


P = 0.764 


5.46(4) X 10-3 


1.53(3) X 10-3 


5.9(2) X lO-'' 


3.0(2) X lO-'' 


2.0(2) X 10-^ 


1.5(2) X 10-^ 


1.2(2) X 10-^ 


f3 = 0.78 


4.28(3) X 10-3 


1.09(2) X 10-3 


3.7(2) X 10-^ 


1.7(1) X 10-"^ 


1.0(2) X 10-^ 


7(1) X 10-^ 


6(1) X 10-5 


/3 = 0.80 


3.39(2) X 10~3 


8.4(2) X lO"'' 


2.7(1) X 10-^ 


1.2(1) X 10-"^ 


6.6(8) X 


4.6(7) X 10-5 


3.6(7) X 10-5 


P = 0.88 


1.441(4) X 10-3 


2.40(2) X 10"^ 


4.87(6) X 10"^ 


1.14(3) X 10-^ 


3.0(1) X 10-6 




3.3(4) X lO-'^ 



TABLE IV: Fitted residual mass values for valence quarks with different extents of the fifth dimension for the Iwasaki gauge action on each 
of our Ls = 8 dynamical ensembles. The valence and sea quark masses were mud = rUg = 0.04. 



Dataset 


Ls = 8 


i^. = 12 


Ls = lQ 


Ls = 20 


i^. = 24 


Ls = 32 


/3 = 2.13 


1.18(1) X 10-2 


4.53(5) X 10-3 


2.25(5) X 10-3 


1.38(4) X 10-3 


9.8(4) X 10" 


-"^ 6.3(3) X lO-'' 


/? = 2.2 


7.19(6) X 10-^ 


2.25(4) X 10-3 


9.0(3) X 10-^ 


4.6(2) X 10-^ 


2.9(2) X 10" 


-4 1.7(2) X 10-4 


/3 = 2.3 


3.90(2) X 10-3 


9.2(1) X 10""^ 


2.65(6) X 10-"^ 


9.5(5) X 10-5 


4.2(4) X 10" 


-5 2.2(4) X 10-5 



TABLE V: Fit parameters to the dependence of mres on for the Iwasaki and DBW2 gauge 
actions. We also give the lattice spacing obtained from the static potential S^] on the dynamical 
background for my^^ = nig = 0.04 without chiral extrapolation. This differs from the lattice spacing 
defined in the chiral limit presented in and we only provide the lattice spacings here for 
illustration. 



Action 


/3 


Cl 


Ac 


C2 


a-i(GeV) 


Iwasaki 


2.13 


0.344(6) 


0.191(2) 


0.0198(3) 


1.62(2) 


Iwasaki 


2.2 


0.301(4) 


0.219(2) 


0.0054(1) 


1.89(2) 


Iwasaki 


2.3 


0.262(3) 


0.268(2) 


0.0006(1) 


2.26(7) 


DBW2 


0.72 


0.342(10) 0.201(4) 


0.0220(4) 


1.40(4) 


DBW2 


0.764 


0.298(6) 


0.252(3) 


0.0040(1) 


1.74(3) 


DBW2 


0.78 


0.267(5) 


0.264(2) 


0.00188(7) 


1.82(2) 


DBW2 


0.80 


0.216(5) 


0.265(3) 


0.00120(5) 


1.98(4) 


DBW2 


0.88 


0.171(7) 


0.338(5) 


1.0(5) X 10-^ 
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0, <\|/'yv|/>^0 X >0 

FIG. 1: Schematic diagram of the quenched Aoki phase. A pionic condensate and corresponding 
non-zero density of near zero modes is developed throughout most of the negative mass region. 
In the coloured sectors, a non-zero mobility edge is developed and the contours could equally well 
represent either decreasing low mode density, decreasing pionic condensate or increasing Ac as one 
moves towards the continuum limit at = 0. 
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U.3 




FIG. 2: Spectral flow of DBW2. The horizontal axis is the valence domain wall height, < M5 < 
3.0 and the vertical axis is the eigenvalue of Hw, —0.3 < < 0.3. From left, f3 = 0.72, ,/3 = 
0.764, /3 = 0.78 (first row) and /? = 0.8,13 = 0.88 (second row). 
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FIG. 3: Spectral flow of Iwasaki /? = 2.13, /3 = 2.2, and /3 = 2.3 (from left to right). Axes are the 
same as fig [2j 
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0.2 
lambda 



FIG. 4: Spectral density obtained on a subset of our data with no error analysis, and fixed bin width 
0.02. Thus, where the density is low the error could be large and the data should be considered 
only as a qualitative indication of the nature of our ensembles. 
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logcp 




FIG. 5: Various measures can be used to define the degree of localization. Inverse participation 
ratio has been a popular order parameter for delocalization transitions and, while fine for that 
purpose, counts occupied sites without paying attention to mode shape. This one dimensional 
cartoon highlights the advantages of the more robust measure employed in this paper. Given the 
model of exponential localization, our tactic is to find a localization exponent that forms a tight 
bound for the eigenvector. This measure will show the influence of satellite peaks in a robust 
fashion by eliminating all forms of dilution by volume average. Other useful measures that take 
into account the shape would include appropriately weighted moments of the density function. 
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DBW2 0.764 



0.1 0.2 0.3 0.4 

Lambda 

FIG. 6: Microscopic view of the mobility edge based on a scatter plot of maximal localization 
lengths of individual low eigenmodes of Hw- We overlay the Hj- mobility edge Ac from fits to 
the model for mres, Eq[35]and find good agreement with our diverging localisation length for the 
DBW2 gauge action at /? = 0.764. 
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Iwasaki2.13 



Lambda 



FIG. 7: Microscopic view of the mobility edge based on a scatter plot of maximal localization 
lengths of individual low eigenmodes of Hw- We overlay the Hj- mobility edge Ac from fits to 
the model for mres, Eq[35]and find good agreement with our diverging localisation length for the 
Iwasaki gauge action at P = 2.13. 
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FIG. 8: We display an iso-surface through the eigenvector density of a typical near-zero chiral 
mode of the 5-dimensional operator //dwf, bound to a single wall. Here s is the vertical direction, 
t the left-right axis, and z runs into the page. We show (x, y) = (3, 3) showing a localized incursion 
of the eigenmode into the fifth dimension associated with a low mode of Hw. This mechanism is 
responsible for one form of the contributions to chiral symmetry breaking in DWF. 
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FIG. 9: We display an iso-surface through the eigenvector density of the same chiral mode displayed 
in FigO Here s is the vertical direction, t the left-right axis, and z runs into the page. We show 
(x,y) = (11, 13) showing a uniformly well stuck region of this mode. Compared with the previous 
figure, this degree of binding is the more common case for this typical eigenvector, and one has to 
search spatially to find any localised incursion. 
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FIG. 10: Normalized overlap of different s-slices of a chiral 5-dimensional eigenmode of H^wp 
with the first 256 4-dimensional eigenmodes of Hyy. We take the valence operator with Lg = 16, 
and demonstrate numerically that only a few low-lying 4-dimensional eigenmodes dominate the 
large Lg contributions to chirality mixing operators such as rrircs by plotting the decomposition 
of the lowest 5-dimensional mode |aoj(s)p as a function of fifth coordinate s, and 4-dimensional 
mode number j. 
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FIG. 11: Residual mass fits for valence quarks with m = 0.04 and Lg = 8, 12, 16,24,32 on the 
DBW2 P = 0.72 ensemble with sea quark masses m^^ = rUg = 0.04, and for the sea quark Ls = 8. 



49 



0.001 



Pi 



0.0001 



le-05 



^ ^ ^ ^ 



3: -- 



I 



-X- 



^ X T ^ 
4 1^^11 



I 



10 
T 



15 



20 



FIG. 12: Residual mass fits for valence quarks witli m = 0.04 and Lg = 8, 12, 16, 24, 32 on the 
DBW2 P = 0.764 ensemble with sea quark masses rriud = rUg = 0.04 and for the sea quark Lg = 8. 
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FIG. 13: Residual mass fits for valence quarks with m = 0.04 and Lg = 8, 12, 16, 24, 32 on the 
DBW2 P = 0.78 ensemble with sea quark masses rriud = = 0.04 and for the sea quark Lg = 8. 
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FIG. 14: Residual mass fits for valence quarks witli m = 0.04 and Lg = 8, 12, 16, 24, 32 on the 
DBW2 P = 0.80 ensemble with sea quark masses rriud = = 0.04 and for the sea quark Lg = 8. 



52 



o.ou 



0.001 



0.0001 



le-05 



le-06 



+ — + + + + + + \ + + + + f 



— + — + t — + — + — + + + — + — f 



T 



10 
T 



15 



20 



FIG. 15: Residual mass fits for valence quarks with m = 0.04 and Lg = 8, 12, 16,24,32 on the 
DBW2 P = 0.88 ensemble with sea quark masses niud = = 0.04 and for the sea quark = 8. 
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FIG. 16: Residual mass fits for valence quarks witli m = 0.04 and Lg = 8, 12, 16, 24, 32 on the 
Iwasaki /? = 2.13 ensemble with sea quark masses m^rf = = 0.04 and for the sea quark Lg = 8. 
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FIG. 17: Residual mass fits for valence quarks witli m = 0.04 and Lg = 8, 12, 16, 24, 32 on the 
Iwasaki /? = 2.2 ensemble with sea quark masses rriy^d = nig = 0.04 and for the sea quark Lg = 8. 
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FIG. 18: Residual mass fits for valence quarks with m = 0.04 and Ls = 8, 12, 16, 24, 32 on the 
Iwasaki (3 = 2.3 ensemble with sea quark masses m^,^ = nis = 0.04 and for the sea quark Lg = 8. 
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FIG. 19: Valence rUj-Qs^Lg) dependence of partially quenched DWF with the DBW2 gauge action 
as a function of L^, using Lg = 8 and the indicated /3 values on a 16'^ x 32 volume for the ensembles. 
The sea and valence quark masses are 0.04 throughout. 
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FIG. 20: Valence mres(-^s) dependence of partially quenched DWF with the Iwasaki gauge action 
as a function of L^, using Lg = 8 and the indicated /? values on a 16'^ x 32 volume for the ensembles. 
The sea and valence quark masses are 0.04 throughout. 
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FIG. 21: Dependence of the plaquette as a function of the dynamical L^, obtained from low 
statistics runs. The dependence of this observable shows the expected C + e~"^' for an exponen- 
tial approach to = oo, and thus demonstrates nicely the connection between finite Lg DWF 
simulations and those with exactly chiral fermions. 
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FIG. 22: The Monte Carlo trajectory history for global topological charge with the DBW2 gauge 
action at couplings (3 = 0.72,0.764,0.78,0.8. The (3 = 0.88 history is offset by 3000 and displays 
a clear loss of topology tunneling. With a reasonable computational effort, a simulation with the 
current algorithm and these parameters would be practically non-ergodic. 
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FIG. 23: The Monte Carlo trajectory history for gfobal topological charge with the Iwasaki gauge 
action at couplings (3 = 2.13,2.2,2.3. The tunneling rate remains acceptable and leaves headroom 
for taking a continuum limit with acceptable topological sampling without any further algorithmic 
improvement. 
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